gezeitenreibung_final.ipynb möglichst nicht ändern, sondern nur die Datei gezeitenreibung.ipynb bearbeiten. Nur fertige komplette Sachen in gezeitenreibung_final.ipynb kopieren.
--
Zellen von gezeitenreibung.ipynb in gezeitenreibung_final.ipynb kopieren. $ \newline $ Baryzentrum in barycenter umbenenen, das wird aber im Augenblick zu verstreut verwendet.
$ \def\dv#1{ \frac{\mathrm{d}}{\mathrm{d} #1} } $ $ \def\dd#1{ \ \mathrm{d} #1 } $ $ \def\Vec#1{ \overrightarrow{#1} } $ $ \def\VecT#1#2{ \begin{pmatrix} #1 \\ #2 \end{pmatrix} } $ $ \def\VecD#1#2#3{ \begin{pmatrix} #1 \\ #2 \\ #3 \end{pmatrix} } $ $ \def\VecE#1{ \hat{e}_{#1} } $ $\providecommand{\e}[1]{\ensuremath{\cdot 10^{#1}}}$ $\providecommand{\fehlt}{\textcolor{red}{\textbf{Fehlt!\dots}}}$ $\providecommand{\todo}{\textcolor{red}{\textbf{\huge{ToDo}}}}$
mErde, mOzean, RErde, mMond, rMondBahn, TMondBahn, G, TErdRotation, tau = 5.9721986*10**24, 0.0014*10**24, 6.3675*10**6, 7.3459*10**22, 3.836*10**8, 27.32166140*24*3600, 6.67430*10**(-11), 86164.100, 0.0021
import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import pi
import scipy
from scipy.integrate import solve_ivp
from scipy.integrate._ivp.ivp import OdeResult # :) https://github.com/scipy/scipy/blob/main/scipy/integrate/_ivp/ivp.py
from matplotlib.animation import FuncAnimation
from IPython.display import Image #to display animations
from matplotlib import cm
# Beautiful plots
plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = "serif"
plt.rcParams['font.serif'] = "Computer Modern Roman"
# Useful StackPosts
# https://stackoverflow.com/questions/59634279/solve-ivp-error-required-step-size-is-less-than-spacing-between-numbers
# Executing everything all the time takes too long
fastExecution = True
# NICHT ANFASSEN! WENN, DANN UNBEDINGT IN FINAL AUCH ÄNDERN!
def baryzentrum(m1, m2, r): return r*m2/(m1+m2)
def iv_stable_orbit_2body():
abstand_baryzentrum_erde = baryzentrum(mErde, mMond, rMondBahn)
return [[-abstand_baryzentrum_erde, 0], [rMondBahn-abstand_baryzentrum_erde, 0], [0, -2*pi*abstand_baryzentrum_erde/TMondBahn], [0, 2*pi*(rMondBahn-abstand_baryzentrum_erde)/TMondBahn], [mErde, mMond]]
def eq_motion_2body(t, state, mass):
x_1, y_1, x_2, y_2, vx_1, vy_1, vx_2, vy_2 = state
dist_em = ((x_2 - x_1)**2 + (y_2 - y_1)**2 )**0.5
ax_1 = G * mass[1] / (dist_em**3) * (x_2 - x_1)
ay_1 = G * mass[1] / (dist_em**3) * (y_2 - y_1)
ax_2 = G * mass[0] / (dist_em**3) * (x_1 - x_2)
ay_2 = G * mass[0] / (dist_em**3) * (y_1 - y_2)
return np.array([vx_1, vy_1, vx_2, vy_2, ax_1, ay_1, ax_2, ay_2])
def two_body_problem(pos_body_1: list, pos_body_2: list, vel_body_1: list, vel_body_2: list, mass: list, t_max: float, steps=1, method='RK45', t_start=0):
solution = solve_ivp(eq_motion_2body, [t_start, t_max], [*pos_body_1, *pos_body_2, *vel_body_1, *vel_body_2], args=(mass,), atol=1e-4*rMondBahn, max_step=TMondBahn/1000, steps=steps, method=method)
x1, y1, x2, y2, v_x1, v_y1, v_x2, v_y2 = solution.y
return [solution.t, x1, y1, x2, y2]
def tide_acceleration(x_M, y_M, x_i, y_i, phi_i, ax_E, ay_E, omega_i, m_M):
part_a, part_b = G*m_M* ((y_M - y_i)*np.cos(phi_i) - (x_M - x_i)*np.sin(phi_i)) / ((x_M - x_i)**2 + (y_M - y_i)**2)**(3/2), (RErde*omega_i**2 - (ax_E*np.cos(phi_i) + ay_E*np.sin(phi_i)))
ax_i, ay_i = -part_a*np.sin(phi_i) - part_b*np.cos(phi_i), part_a*np.cos(phi_i) - part_b*np.sin(phi_i)
return (-(ax_i - ax_E)*np.sin(phi_i) + (ay_i - ay_E)*np.cos(phi_i))/RErde
def eq_motion_4body(t, state, mass, calculate_tides=True):
x_E, y_E, x_M, y_M, phi_1, phi_2, vx_E, vy_E, vx_M, vy_M, omega_1, omega_2 = state
m_E, m_M, m_1, m_2 = mass
x_1, y_1, x_2, y_2 = x_E + RErde*np.cos(phi_1), y_E + RErde*np.sin(phi_1), x_E + RErde*np.cos(phi_2), y_E + RErde*np.sin(phi_2)
mu_x = m_E + m_1*np.cos(phi_1)**2 + m_2*np.cos(phi_2)**2
mu_y = m_E + m_1*np.sin(phi_1)**2 + m_2*np.sin(phi_2)**2
a = (G*m_M*m_E)/((x_M - x_E)**2 + (y_M - y_E)**2)**(3/2)
b = m_1*np.sin(phi_1)*np.cos(phi_1) + m_2*np.sin(phi_2)*np.cos(phi_2)
c1 = m_1*RErde*omega_1**2 + (G*m_M*m_1*((x_M - x_1)*np.cos(phi_1) + (y_M - y_1)*np.sin(phi_1)))/((x_M - x_1)**2 + (y_M - y_1)**2)**(3/2)
c2 = m_2*RErde*omega_2**2 + (G*m_M*m_2*((x_M - x_2)*np.cos(phi_2) + (y_M - y_2)*np.sin(phi_2)))/((x_M - x_2)**2 + (y_M - y_2)**2)**(3/2)
ax_E = 1/(mu_x - b**2/mu_y) * (a*(x_M - x_E) -b/mu_y*(a*(y_M - y_E) + c1*np.sin(phi_1) + c2*np.sin(phi_2)) + c1*np.cos(phi_1) + c2*np.cos(phi_2))
ay_E = 1/mu_y*(a*(y_M - y_E) - b*ax_E + c1*np.sin(phi_1) + c2*np.sin(phi_2))
ax_M = G*((m_E*(x_E - x_M))/((x_E - x_M)**2 + (y_E - y_M)**2)**(3/2) + (m_1*(x_1 - x_M))/((x_1 - x_M)**2 + (y_1 - y_M)**2)**(3/2) + (m_2*(x_2 - x_M))/((x_2 - x_M)**2 + (y_2 - y_M)**2)**(3/2))
ay_M = G*((m_E*(y_E - y_M))/((x_E - x_M)**2 + (y_E - y_M)**2)**(3/2) + (m_1*(y_1 - y_M))/((x_1 - x_M)**2 + (y_1 - y_M)**2)**(3/2) + (m_2*(y_2 - y_M))/((x_2 - x_M)**2 + (y_2 - y_M)**2)**(3/2))
if not calculate_tides: return [vx_E, vy_E, vx_M, vy_M, omega_1, omega_2, ax_E, ay_E, ax_M, ay_M, 0, 0] # For code reuse with friction
alpha_1, alpha_2 = tide_acceleration(x_M, y_M, x_1, y_1, phi_1, ax_E, ay_E, omega_1, m_M), tide_acceleration(x_M, y_M, x_2, y_2, phi_2, ax_E, ay_E, omega_2, m_M)
return[vx_E, vy_E, vx_M, vy_M, omega_1, omega_2, ax_E, ay_E, ax_M, ay_M, alpha_1, alpha_2]
def center_of_mass(mass, position_x , position_y):
return [np.sum(mass*position_x)/np.sum(mass), np.sum(mass*position_y)/np.sum(mass)]
def iv_stable_orbit_4body():
mass = [mErde, mMond, mOzean/2, mOzean/2] # Each of the tides gets half of the ocean's mass
abstand_baryzentrum_erde = center_of_mass(mass, np.array([0, rMondBahn, -RErde, RErde]), np.array([0, 0, 0, 0]))[0]
x0_Erde = [-abstand_baryzentrum_erde, 0]
x0_Mond = [rMondBahn - abstand_baryzentrum_erde, 0]
vy0_Erde = 2*pi*abstand_baryzentrum_erde/TMondBahn
vy0_Mond = 2*pi*(rMondBahn-abstand_baryzentrum_erde)/TMondBahn
v0_Erde = [0, -vy0_Erde]
v0_Mond = [0, vy0_Mond]
phi_1, phi_2 = 0, pi
omega_1, omega_2 = 2*pi/TMondBahn, 2*pi/TMondBahn
return [x0_Erde, x0_Mond, [phi_1, phi_2], v0_Erde, v0_Mond, [omega_1, omega_2], mass]
def four_body_problem(pos_body_1: list, pos_body_2: list, phi_i: list, vel_body_1: list, vel_body_2: list, omega_i: list, mass: list, t_max: float, fun=eq_motion_4body, rtol=1e-12, atol=1e-12, t_start=0):
solution = solve_ivp(fun=fun, t_span=[t_start, t_max], y0=[*pos_body_1, *pos_body_2, *phi_i, *vel_body_1, *vel_body_2, *omega_i], args=(mass,), rtol=rtol, atol=atol)
x_E, y_E, x_M, y_M, phi_1, phi_2, vx_E, vy_E, vx_M, vy_M, omega_1, omega_2 = solution.y
return [solution.t, x_E, y_E, x_M, y_M, phi_1, phi_2, omega_1, omega_2]
# Calculate the solution of the two body problem
t, x_1, y_1, x_2, y_2 = two_body_problem(*iv_stable_orbit_2body(), 3*TMondBahn)
/home/finn/.local/lib/python3.10/site-packages/scipy/integrate/_ivp/common.py:39: UserWarning: The following arguments have no effect for a chosen solver: `steps`.
warn("The following arguments have no effect for a chosen solver: {}."
Eine exakte analytische Lösung des Zweikörperproblems ist mit genügend Annahmen möglich. Dazu wird im Schwerpunktssytem gearbeitet. Aus den Bewegungsgliechungen der Körper die nur mit der Schwerkraft wechselwirken, kann man folgern, dass der Schwerpunkt des Systems eine gleichförmige Bewegung ausführt. Mit unseren Anfangsbedingungen bewegt dieser sich aber nicht.
Die beiden Körper können sich um ihren Schwerpunkt je nach Anfangsbedingungen auf verschiedenen Bahnen bewegen. Dies lässt sich durch die Exzentrizität beschreiben. Mögliche Bahnen sind Kreise, Ellipsen, Parabeln und Hyperbeln. Wir haben die Anfangsbedingungen so gewählt, dass die Bahn eine Kreisbahn ist. $\newline$ Unsere Ergebnise stimmen also mit der exakten Lösung des Zweikörperproblems überein, da wir zwei geschlossene Kreisbahnen um den gemeinsamen Schwerpunkt (mit unseren Anfangsbedingungen der Ursprung) erhalten. Diese Kreisbahnen lassen sich mit einfachen trigonometrischen Funktionen beschreiben. \begin{align} r_E(t) = \VecT{x_E(t)}{y_E(t)} = \VecT{r_{E_0}\cos(\omega_E t)}{r_{E_0}\sin(\omega_E t)} \\ r_M(t) = \VecT{x_M(t)}{y_M(t)} = \VecT{r_{M_0}\cos(\omega_M t)}{r_{M_0}\sin(\omega_M t)} \end{align} Diese lassen sich ebenfalls leicht plotten.
Unser numerischer Ansatz ist natürlich in der Lage alle möglichen Bahnen zu berechnen, also auch Ellipsen und Hyperbeln, es ließe sich also auch ein Komet oder ein Sling-Shot-Manöver berechnen. Mit den geeigneten Anfangsbedingungen ließen sich diese auch analytisch berechnen.
Eine weitere Möglichkeit unsere Ergebnisse zu überprüfen ist das dritte Kepplersche Gesetz, es besagt, dass die Umlaufzeit $T$ des Zweikörpersystem beträgt: $$T^2 = \frac{4\pi^2 \cdot a^3}{G\cdot M} $$ Mit $G$ als Gravitationskonstante, $M = m_{Erde}+m_{Mond}$ und $a$ als großer Halbachse welche in unserem genäherten Fall der Kreisbahn $a=r_{MondBahn}$ entspricht.
r_E0, r_M0 = iv_stable_orbit_2body()[0][0], iv_stable_orbit_2body()[1][0]
time_points = np.linspace(0, 3*TMondBahn, 1000)
x_E, y_E = np.vectorize(lambda t: r_E0 * np.cos(2*pi/TMondBahn * t))(time_points), np.vectorize(lambda t: r_E0 * np.sin(2*pi/TMondBahn * t))(time_points)
x_M, y_M = np.vectorize(lambda t: r_M0 * np.cos(2*pi/TMondBahn * t))(time_points), np.vectorize(lambda t: r_M0 * np.sin(2*pi/TMondBahn * t))(time_points)
# Plot the analytical solution with our boundary conditions
plt.rc('axes', axisbelow=True)
plt.figure(figsize=(5, 5))
plt.title("Analytische Bahnen des 2-Körper-Systems")
plt.grid('minor', 'minor', linestyle='-', linewidth=0.2)
plt.grid('major', 'major', linestyle='-', linewidth=0.8)
plt.minorticks_on()
plt.xlabel("x")
plt.ylabel("y")
plt.plot(x_E, y_E, 'b', linewidth=1, label='Erde')
plt.plot(x_M, y_M, 'r', linewidth=1, label='Mond')
plt.legend()
plt.show()
print(f'3. Keplersches Gesetze T: {TMondBahn} = {(4*pi**2*rMondBahn**3/(G*(mErde+mMond)))**(1/2)}')
3. Keplersches Gesetze T: 2360591.54496 = 2350028.5085322363
Im Folgenden wollen wir das 4-Körper-Problem um die Reibung des Ozeans (der Flutberge) mit der Erde erweitern. Dies verlangsamt die intrinsische Erdrotation und sorgt für länger werdende Tage und eine größere Umlaufbahn des Mondes. Um dies zu beschreiben müssen wir mehrere zusätzliche Komponenten herleiten und mit ihnen unser Modell erweitern.
Die Reibung beschreiben wir zunächst mit der Reibungskonstante $ k = 2 \cdot 10^{-12}\,\frac{1}{\mathrm m}$.
Wir wollen als erstes die Reibungskraft zwischen der Erde und den Flutbergen bestimmen. Wir wählen dafür folgenden Ansatz: $ \vec F_{\mathrm R,i}(\vec v_i)=-km_i|\vec v_i| \vec v_i $. Hierbei ist $\vec v_i$ die Geschwindigkeit des $i$-ten Flutbergs relativ zur rotierenden Erdoberfläche und $k$ eine effektive Reibungskonstante.
Die Geschwindigkeit des Flutbergs bestimmen wir über seine Winkelgeschwindikeit relativ zur Erdoberfläche $(\omega_E - \omega_i)$, mit der er sich um die z-Achse dreht. Hierbei ist $\omega_E $ die Winkelgeschwindigkeit der Eigenrotation der Erde und $\omega_i$ die Winkelgeschwindigkeit des Flutbergs.
Die Geschwindigkeit ist allgemein definiert als $\vec v_i = \vec \omega_i \times \vec r_i$: \begin{equation*} \vec v_i = (\omega_E - \omega_i) \cdot \VecD{0}{0}{1} \times r_{E} \cdot \vec e_{r, i} = \VecD{0}{0}{\omega_E - \omega_i} \times \VecD{r_{E} \cdot \cos(\phi_i)}{r_{E} \cdot \sin(\phi_i)}{0} = \VecD{-r_{E} \cdot \sin(\phi_i) \cdot (\omega_E - \omega_i)}{r_{E} \cdot \cos(\phi_i) \cdot (\omega_E - \omega_i)}{0} \end{equation*} Damit berechnen wir den Betrag: \begin{equation*} |\vec v_i| = \sqrt{r_{E}^2 \cdot (\omega_E - \omega_i)^2 \cdot \sin^2(\phi_i) + r_{E}^2 \cdot (\omega_E - \omega_i)^2 \cdot \cos^2(\phi_i)} = r_{E} \cdot |\omega_E - \omega_i| \end{equation*} Also ist die Reibungskraft: \begin{equation*} \vec F_{\mathrm R,i}(\vec v_i) = -km_i |\vec v_i| \vec v_i = -km_i r_{E} \cdot |\omega_E - \omega_i| \cdot \VecD{-r_{E} \cdot \sin(\phi_i) \cdot (\omega_E - \omega_i)}{r_{E} \cdot \cos(\phi_i) \cdot (\omega_E - \omega_i)}{0} = -km_i r_{E}^2 \cdot |\omega_E - \omega_i|(\omega_E - \omega_i) \cdot \VecD{- \sin(\phi_i)}{\cos(\phi_i)}{0} \end{equation*} Durch teilen durch die Masse des i.ten Flutbergs erhalten wir die Beschleuningung, die wir auf die Differentialgleichung des i.ten Flutbergs aus dem 4-Körper-Problem addieren.
Als nächstes wollen wir die Differentialgleichung für die intrinsische Rotation der Erde herleiten. Es gilt die Drehimpulsrelation $L = I \cdot \omega$. Wir nehmen die Erde als Kugel mit homogener Massenverteilung an, damit hat sie ein Trägheitsmoment von $I = \frac{2}{5} m_E r_E^2$. Für Drehimpuls und Drehmoment folgt: \begin{gather*} \vec L = \frac{2}{5} m_E r_E^2 \vec \omega_E \\ \vec M = \frac{\mathrm d \vec L}{\mathrm d t} = \frac{2}{5} m_E r_E^2 \vec{\dot \omega_E} \end{gather*} Das Drehmoment welches durch die Reibungskräfte entsteht ist: \begin{equation*} \vec M_{\mathrm R} = \sum_{i=1}^2 \vec r_i \times \vec F_{\mathrm R,i} \end{equation*} Dies berechnet sich für zwei Flutberge zu: \begin{equation*} \vec M_{\mathrm R} = \VecD{0}{0}{-k r_E^3 \left(m_1 \cdot |\omega_E - \omega_1|(\omega_E - \omega_1) + m_2 \cdot |\omega_E - \omega_2|(\omega_E - \omega_2) \right)} \end{equation*} Setzen wir die z-Komponenten der Drehmomente gleich (nur diese sind $\neq 0$), erhalten wir folgende Bewegungsgleichung für die intrinsische Rotation der Erde: \begin{gather*} \frac{2}{5} m_E r_E^2 \dot \omega_E = -k r_E^3 \left(m_1 \cdot |\omega_E - \omega_1|(\omega_E - \omega_1) + m_2 \cdot |\omega_E - \omega_2|(\omega_E - \omega_2) \right) \\ \dot \omega_E = -\frac{5 k r_E}{2 m_E} \left(m_1 \cdot |\omega_E - \omega_1|(\omega_E - \omega_1) + m_2 \cdot |\omega_E - \omega_2|(\omega_E - \omega_2) \right) \end{gather*}
Um dies im Code umzusetzen verwenden wir, wenn es geht, die schon implementierten Funktionen. Daher wird hier der selbe Ansatz wie zuvor gewählt, nur die Beschleunigung auf die Flutberge wird wie oben erklärt erweitert, und die zusätzliche Differentialgleichung wird hinzugefügt.
def iv_stable_orbit_4body_friction():
'''Initial conditions for the Earth-Moon-Tides system in a stable orbit'''
x0_Erde, x0_Mond, phi_i, v0_Erde, v0_Mond, omega_i, mass = iv_stable_orbit_4body()
phi0_E = 0 # The starting point of the earth's rotation is unconsequentially defined as 0
omega0_E = 2*pi / TErdRotation # Angular velocity of the earth
return [x0_Erde, x0_Mond, phi_i, phi0_E, v0_Erde, v0_Mond, omega_i, omega0_E, mass]
def tides_acceleration_friction(x_M, y_M, x_i, y_i, phi_i, ax_E, ay_E, omega_i, omega_E, m_M, k):
'''Calculates the angular acceleration of a tide, adding the frictional force caused by the earth's rotation'''
ax_friction = k * RErde**2 * np.abs(omega_i - omega_E)* (omega_i - omega_E) * np.sin(phi_i)
ay_friction = -k * RErde**2 * np.abs(omega_i - omega_E)* (omega_i - omega_E) * np.cos(phi_i)
part_a = G*m_M* ((y_M - y_i)*np.cos(phi_i) - (x_M - x_i)*np.sin(phi_i)) / ((x_M - x_i)**2 + (y_M - y_i)**2)**(3/2)
part_b = (RErde*omega_i**2 - (ax_E*np.cos(phi_i) + ay_E*np.sin(phi_i)))
ax_i = -part_a*np.sin(phi_i) - part_b*np.cos(phi_i) + ax_friction
ay_i = part_a*np.cos(phi_i) - part_b*np.sin(phi_i) + ay_friction
alpha_i = (-(ax_i - ax_E)*np.sin(phi_i) + (ay_i - ay_E)*np.cos(phi_i))/RErde
return alpha_i
def eq_motion_4body_friction(t, state, mass, k):
x_E, y_E, x_M, y_M, phi_1, phi_2, phi_E, vx_E, vy_E, vx_M, vy_M, omega_1, omega_2, omega_E = state
x_1, y_1, x_2, y_2 = x_E + RErde*np.cos(phi_1), y_E + RErde*np.sin(phi_1), x_E + RErde*np.cos(phi_2), y_E + RErde*np.sin(phi_2)
m_E, m_M, m_1, m_2 = mass
# Reuse the equations of motion from the 4-body problem for the earth and the moon:
_, _, _, _, _, _, ax_E, ay_E, ax_M, ay_M, _, _ = eq_motion_4body(t, np.array([x_E, y_E, x_M, y_M, phi_1, phi_2, vx_E, vy_E, vx_M, vy_M, omega_1, omega_2]), mass, calculate_tides=False)
# Tides angular acceleration with friction:
alpha_1 = tides_acceleration_friction(x_M, y_M, x_1, y_1, phi_1, ax_E, ay_E, omega_1, omega_E, m_M, k)
alpha_2 = tides_acceleration_friction(x_M, y_M, x_2, y_2, phi_2, ax_E, ay_E, omega_2, omega_E, m_M, k)
# The earth's angular acceleration:
alpha_E = (5*k*RErde)/(2*m_E) * (m_1*np.abs(omega_1 - omega_E)*(omega_1 - omega_E) + m_2*np.abs(omega_2 - omega_E)*(omega_2 - omega_E))
return[vx_E, vy_E, vx_M, vy_M, omega_1, omega_2, omega_E, ax_E, ay_E, ax_M, ay_M, alpha_1, alpha_2, alpha_E]
def four_body_problem_friction(pos_body_1: list, pos_body_2: list, phi_i: list, phi_E: float, vel_body_1: list, vel_body_2: list, omega_i: list, omega_E: float, mass: list, k: float, t_max: float, rtol=1e-6, atol=1e-6, t_start=0):
solution = solve_ivp(fun=eq_motion_4body_friction, t_span=[t_start, t_max], y0=[*pos_body_1, *pos_body_2, *phi_i, phi_E, *vel_body_1, *vel_body_2, *omega_i, omega_E], args=(mass, k), rtol=rtol, atol=atol)
x_E, y_E, x_M, y_M, phi_1, phi_2, phi_E, vx_E, vy_E, vx_M, vy_M, omega_1, omega_2, omega_E = solution.y
return [solution.t, x_E, y_E, x_M, y_M, phi_1, phi_2, phi_E, omega_1, omega_2, omega_E]
if not fastExecution:
# 4 body without friction:
t, x_E, y_E, x_M, y_M, phi_1, phi_2, _, _ = four_body_problem(*iv_stable_orbit_4body(), 0.75*TMondBahn, rtol=1e-12, atol=1e-12) # 0.3s
x_1, y_1, x_2, y_2 = x_E + RErde*np.cos(phi_1), y_E + RErde*np.sin(phi_1), x_E + RErde*np.cos(phi_2), y_E + RErde*np.sin(phi_2)
t_l, x_E_l, y_E_l, x_M_l, y_M_l, phi_1_l, phi_2_l, omega_1_l, omega_2_l = four_body_problem(*iv_stable_orbit_4body(), 100*365*24*3600, rtol=1e-9, atol=1e-9)
x_1_l, y_1_l, x_2_l, y_2_l = x_E_l + RErde*np.cos(phi_1_l), y_E_l + RErde*np.sin(phi_1_l), x_E_l + RErde*np.cos(phi_2_l), y_E_l + RErde*np.sin(phi_2_l)
# 4 body with friction short time: 0.3s
t_f, x_E_f, y_E_f, x_M_f, y_M_f, phi_1_f, phi_2_f, phi_E_f, omega_1_f, omega_2_f, omega_E_f = four_body_problem_friction(*iv_stable_orbit_4body_friction(), k=2e-12, t_max=0.75*TMondBahn, rtol=1e-12, atol=1e-12)
x_1_f, y_1_f, x_2_f, y_2_f = x_E_f + RErde*np.cos(phi_1_f), y_E_f + RErde*np.sin(phi_1_f), x_E_f + RErde*np.cos(phi_2_f), y_E_f + RErde*np.sin(phi_2_f)
# 4 body with friction for 100 years: 11min
t_fl, x_E_fl, y_E_fl, x_M_fl, y_M_fl, phi_1_fl, phi_2_fl, phi_E_fl, omega_1_fl, omega_2_fl, omega_E_fl = four_body_problem_friction(*iv_stable_orbit_4body_friction(), k=2e-12, t_max=100*365*24*3600, rtol=1e-12, atol=1e-12)
x_1_fl, y_1_fl, x_2_fl, y_2_fl = x_E_fl + RErde*np.cos(phi_1_fl), y_E_fl + RErde*np.sin(phi_1_fl), x_E_fl + RErde*np.cos(phi_2_fl), y_E_fl + RErde*np.sin(phi_2_fl)
if not fastExecution:
# Plot the solution of the four body problem as a phase space diagram
fig, axes = plt.subplots(nrows=5, ncols=4, figsize=(12, 12)) # Create a figure and a set of subplots
fig.suptitle("Phasenraumdiagramme des Erde-Mond-Flutberg-Systems: ohne Reibung und mit Reibung")
def init_phase_space(ax, a, b, title, xlabel, ylabel, color):
ax.plot(a, b, color=color)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(title)
init_phase_space(axes[0,0], t, x_E, "x-Position Erde", "$t$ / s", "$x_1$ / m", "blue")
init_phase_space(axes[0,1], t, x_M, "x-Position Mond", "$t$ / s", "$x_2$ / m", "red")
init_phase_space(axes[1,0], t, y_E, "y-Position Erde", "$t$ / s", "$y_1$ / m", "blue")
init_phase_space(axes[1,1], t, y_M, "y-Position Mond", "$t$ / s", "$y_2$ / m", "red")
init_phase_space(axes[2,0], x_E, y_E, "Rotation Erde", "$x_E$ / m", "$y_E$ / m", "blue")
init_phase_space(axes[2,1], x_M, y_M, "Rotation Mond", "$x_M$ / m", "$y_M$ / m", "red")
init_phase_space(axes[3,0], t, phi_1, "Winkel Flutberg 1", "$t$ / s", "$\phi_1$ / rad", "blue")
init_phase_space(axes[3,1], t, phi_2, "Winkel Flutberg 2", "$t$ / s", "$\phi_2$ / rad", "red")
init_phase_space(axes[4,0], x_1, y_1, "Rotation Flutberg 1", "$x_1$ / m", "$x_2$ / m", "blue")
init_phase_space(axes[4,1], x_2, y_2, "Rotation Flutberg 2", "$y_1$ / m", "$y_2$ / m", "red")
init_phase_space(axes[0,2], t_f, x_E_f, "x-Position Erde mit Reibung", "$t$ / s", "$x_1$ / m", "blue")
init_phase_space(axes[0,3], t_f, x_M_f, "x-Position Mond mit Reibung", "$t$ / s", "$x_2$ / m", "red")
init_phase_space(axes[1,2], t_f, y_E_f, "y-Position Erde mit Reibung", "$t$ / s", "$y_1$ / m", "blue")
init_phase_space(axes[1,3], t_f, y_M_f, "y-Position Mond mit Reibung", "$t$ / s", "$y_2$ / m", "red")
init_phase_space(axes[2,2], x_E_f, y_E_f, "Rotation Erde mit Reibung", "$x_E$ / m", "$y_E$ / m", "blue")
init_phase_space(axes[2,3], x_M_f, y_M_f, "Rotation Mond mit Reibung", "$x_M$ / m", "$y_M$ / m", "red")
init_phase_space(axes[3,2], t_f, phi_1_f, "Winkel Flutberg 1 mit Reibung", "$t$ / s", "$\phi_1$ / rad", "blue")
init_phase_space(axes[3,3], t_f, phi_2_f, "Winkel Flutberg 2 mit Reibung", "$t$ / s", "$\phi_2$ / rad", "red")
init_phase_space(axes[4,2], x_1_f, y_1_f, "Rotation Flutberg 1 mit Reibung", "$x_1$ / m", "$x_2$ / m", "blue")
init_phase_space(axes[4,3], x_2_f, y_2_f, "Rotation Flutberg 2 mit Reibung", "$y_1$ / m", "$y_2$ / m", "red")
plt.tight_layout()
plt.savefig("4body_friction_comp0.png", dpi=300)
plt.close()
fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(15, 10))
plt.rc('axes', axisbelow=True)
plt.suptitle("Verlgeich 4 Körper Problem mit Reibung / ohne Reibung")
def plot_subplot(ax, lines, title, xlabel, ylabel):
for x, y, label, color in lines:
ax.plot(x, y, color, linewidth=1, label=label)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(title)
ax.grid('minor', 'minor', linestyle='-', linewidth=0.2)
ax.grid('major', 'major', linestyle='-', linewidth=0.8)
ax.minorticks_on()
ax.legend()
# Distance from the origin for earth and moon
r_Earth = np.sqrt(x_E_l**2 + y_E_l**2)
r_Moon = np.sqrt(x_M_l**2 + y_M_l**2)
r_Earth_f = np.sqrt(x_E_fl**2 + y_E_fl**2)
r_Moon_f = np.sqrt(x_M_fl**2 + y_M_fl**2)
plot_subplot(axes[0][0], [(x_E_l*10, y_E_l*10, 'Erde (*10)', 'k'), (x_M_l, y_M_l, 'Mond', 'r'), (x_1_l*10, y_1_l*10, 'Flutberg 1 (*10)', 'b'), (x_2_l*10, y_2_l*10, 'Flutberg 2 (*10)', 'b')], "Bahnen des 4-Körper Problems ohne Reibung", "x / m", "y / m")
plot_subplot(axes[0][1], [(t_l, r_Moon, 'Mond', 'r')], "Abstand des Mondes zum Ursprung ohne Reibung", "t / s", "r / m")
plot_subplot(axes[0][2], [(t_l, r_Earth, 'Erde', 'k')], "Abstand der Erde zum Ursprung ohne Reibung", "t / s", "r / m")
plot_subplot(axes[1][0], [(x_E_fl*10, y_E_fl*10, 'Erde (*10)', 'k'), (x_M_fl, y_M_fl, 'Mond', 'r'), (x_1_fl*10, y_1_fl*10, 'Flutberg 1 (*10)', 'b'), (x_2_fl*10, y_2_fl*10, 'Flutberg 2 (*10)', 'b')], "Bahnen des 4-Körper Problems mit Reibung", "x / m", "y / m")
plot_subplot(axes[1][1], [(t_fl, r_Moon_f, 'Mond', 'r')], "Abstand des Mondes zum Ursprung mit Reibung", "t / s", "r / m")
plot_subplot(axes[1][2], [(t_fl, r_Earth_f, 'Erde', 'k')], "Abstand der Erde zum Ursprung mit Reibung", "t / s", "r / m")
plt.tight_layout()
plt.savefig("4body_friction_comp1.png", dpi=300)
plt.close()
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 5))
plt.rc('axes', axisbelow=True)
plt.suptitle("Verlgeich 4 Körper Problem mit Reibung / ohne Reibung")
plot_subplot(axes[0][0], [(t_l, omega_1_l, '$\omega_1$ Flutberg 1', 'k'), (t_l, omega_2_l, '$\omega_2$ Flutberg 2', 'r'), (t_l, [2*pi/TMondBahn]*np.shape(t_l)[0], 'Startwinkelgeschwindigkeit Flutberg', 'g')], "Winkelgeschwindigkeiten der Flutberge ohne Reibung", "t / s", "$\omega$ / rad/s")
plot_subplot(axes[0][1], [(t_l, [2*pi/TErdRotation]*np.shape(t_l)[0], '$\omega_E$ Winkelgeschwindigkeit Erde', 'k')], "Winkelgeschwindigkeit der Erde ohne Reibung", "t / s", "$\omega$ / rad/s")
plot_subplot(axes[1][0], [(t_fl, omega_1_fl, '$\omega_1$ Flutberg 1', 'k'), (t_fl, omega_2_fl, '$\omega_2$ Flutberg 2', 'r'), (t_fl, [2*pi/TMondBahn]*np.shape(t_fl)[0], 'Startwinkelgeschwindigkeit Flutberg', 'g')], "Winkelgeschwindigkeiten der Flutberge mit Reibung", "t / s", "$\omega$ / rad/s")
plot_subplot(axes[1][1], [(t_fl, omega_E_fl, '$\omega_E$ Winkelgeschwindigkeit Erde', 'k'), (t_fl, [2*pi/TErdRotation]*np.shape(t_fl)[0], 'Startwinkelgeschwindigkeit Erde', 'g')], "Winkelgeschwindigkeit der Erde mit Reibung", "t / s", "$\omega$ / rad/s")
plt.tight_layout()
plt.savefig("4body_friction_comp2.png", dpi=300)
plt.close()
# Display image set the size to be 500x500 pixels
display(Image(filename='4body_friction_comp0.png', width=1000, height=2000))
display(Image(filename='4body_friction_comp1.png', width=1000, height=500))
display(Image(filename='4body_friction_comp2.png', width=1000, height=500))
Im Fall ohne Reibung ist $\omega_\mathrm{E}$ natürlich konstant, und im Fall mit Reibung nimmt $\omega_\mathrm{E}$ mit der Zeit ab, das heißt, dass die Tage auf der Erde länger werden.
Hier sieht die Abnahme linear aus, es handelt sich aber um einen exponentiellen Abfall, welcher aufgrund der betrachteten Zeitskala aber linear erscheint.
Zudem ist zu sehen, dass Erde und Mond sich mit Reibung von ihrem gemeinsamen Schwerpunkt wegbewegen, was ohne Reibung nicht passiert.
Das ist eine Folge davon, dass die Erdrotation schneller als die Umlaufbahn des Mondes ist.
Die Flutberge werden aufgrund der Reibung vor die Verbindungslinie zwischen Erde und Mond gedrückt, was die Wirkung der Flutberge auf den Mond ändert.
Der Mond erfährt eine Radialbeschleunigung, wodurch er auf eine weiter entfernte Bahn gelangt.
Vergleicht man die Winkelgeschwindigkeiten der Flutberge, so erkennt man, dass diese in beiden Fällen mit einer schwach modulierten Oszillation starten.
Mit der Zeit schwächt sich dies mit Reibung ab, ohne Reibung bleibt der Effekt zeitlich gleich.
Die Bewegung ist Folge davon, dass die Umlaufbahn nicht perfekt ein Kreis ist (Wie man auch in den Plots vom Abstand zum Uspung erkennen kann), weshalb der Mond die Flutberge zu unterschiedlichen Zeitpunkten unterschiedlich stark beschleunigt.
Mit Reibung wird diese Bewegung abgebremst (Energiedissipation), weshalb die Amplitude der Oszillation abnimmt.
Die Amplitude ist zu Beginn jedoch größer, weil die Flutberge durch die Reibung in Richtung der Erdrotation beschleunigt werden.
Wir stellen noch dar, wie sich der Mond von der Erde entfernt:
#Show the incerase of the moons orbit
if not fastExecution:
#Moon as seen from Earth:
x = x_M_fl - x_E_fl
y = y_M_fl - y_E_fl
fig = plt.figure(figsize=(12,6))
fig.suptitle('Vergrößerung der Mondbahn über 100 Jahre')
ax = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax.set_title('Mond im Erdsystem')
ax.set_xlabel("$x$ / m")
ax.set_ylabel("$y$ / m")
ax.set_aspect('equal')
ax.plot(x, y, color='r', linewidth=0.5, label='Mond')
rect = plt.Rectangle((3.75e8, -1e5), 1.5e7, 1e7, facecolor="black", alpha=0.7, label='Ausschnitt')
ax.add_patch(rect)
ax.legend()
ax2.set_title('Ausschnitt')
ax2.set_xlabel("$x$ / m")
ax2.set_ylabel("$y$ / m")
ax2.set_aspect('equal')
ax2.plot(x, y, 'r', linewidth=0.1, label='Mond')
ax2.scatter(x[0], y[0], c='k', s=10, label='Startpunkt')
ax2.set_xlim([3.8358e8, 3.8391e8])
ax2.set_ylim([-1e5, 1e5])
ax2.legend()
plt.tight_layout()
plt.savefig("4body_friction_Moon.png", dpi=300)
plt.close()
display(Image(filename='4body_friction_Moon.png', width=1000, height=2000))
Für die erstmalige Berechnung der Tageslängenänderung verwenden wir die Reibungskonstante $k = 2 \cdot 10^{-12}\,\frac{1}{\mathrm m}$ und Massen für Flutberge von $ m_1 = m_2 = \frac{1}{2} m_{\mathrm{Ozean}} = \frac{1}{2} \cdot 0.0014 \cdot 10^{24} $
Die Verlangsamung der Erdrotation kann man auch in Realität beobachten. Deshalb vergleichen wir unsere Ergebnisse mit dem Literaturwert für die Tageslängenänderung über 100 Jahre $ \tau = 0.0021 \frac{\mathrm s}{\mathrm{100 a}}$.
Wie im Graphen für die Winkelgeschwindigkeiten der Flutberge zu erkennen ist braucht das System eine gewisse Einschwingzeit. Wir wählen konservativ 200 Jahre als Einschwingzeit (mit Hilfe des Graphen bestimmt). Im folgenden berechnen wir die Tageslängenänderung für 500 Jahre mit 100 Jahr Zeitschritten.
Die Simulation liefert uns die Winkelgeschwindigkeit der Erde $\omega_E$ zu gewissen Zeitpunkten. $ \newline $ Diese kann mit der Formel $\omega_E = \frac{2 \pi}{T} \Rightarrow T = 2 \pi \omega$ in die Tageslänge $T$ umgerechnet werden. Die Differenz der Tageslängen ergibt die Tageslängenänderung $ \tau $.
if not fastExecution:
year_in_seconds = 365*24*3600
t_f5ha, _, _, _, _, _, _, _, _, _, omega_E_f5ha = four_body_problem_friction(*iv_stable_orbit_4body_friction(), k=2e-12, t_max=500*year_in_seconds, rtol=1e-12, atol=1e-12)
angular_velocities = [omega_E_f5ha[i] for i in [np.argmin(np.abs(t_f5ha - j*year_in_seconds)) for j in [0, 100, 200, 300, 400, 500]]]
tau_calc = [2*pi/angular_velocities[i+1] - 2*pi/angular_velocities[i] for i in range(5)]
for i in range (5):
print(f'τ für ({i}ha-{i+1}ha): {tau_calc[i]:.3f} s/ha')
else:
print("τ für (0ha-1ha): 137.090 s/ha\nτ für (1ha-2ha): 137.086 s/ha\nτ für (2ha-3ha): 137.082 s/ha\nτ für (3ha-4ha): 137.078 s/ha\nτ für (4ha-5ha): 137.073 s/ha")
τ für (0ha-1ha): 137.090 s/ha τ für (1ha-2ha): 137.086 s/ha τ für (2ha-3ha): 137.082 s/ha τ für (3ha-4ha): 137.078 s/ha τ für (4ha-5ha): 137.073 s/ha
Das berechnete Ergebnis von ca. $ 137 \frac{\mathrm s}{\mathrm{100 a}}$ weicht stark vom Literaturwert $\tau$ ab.
Es ist deutlich zu erkennen, dass die Tageslängenänderung mit der Zeit abnimmt, da die Reibungskraft mit der Zeit abnimmt. Auf unseren Zeitskalen ist eine Näherung als linear aber noch vertretbar. Die Einschwingzeit scheint für die Berechnung nicht relevant zu sein, da in den ersten 100 Jahren kein deutlich größerer Unterschied des Ergebnisses zu den folgenden Zeitschritten zu erkennen ist.
Da unser berechneter Wert so stark vom Literaturwert abweicht, möchten wir die Parameter unserer Simulation anpassen um den Literaturwert zu erreichen. Dies verwenden wir für die folgenden Fits, um die Laufzeit zu verkürzen in dem wir die Berechnungen nur für die ersten 100 Jahre durchführen.
Die Reibung der Flutberge auf der Erdoberfläche sorgt für eine Verlangsamung der Erdrotation und einem größeren Orbit für den Mond. Wir haben für diese Berechnung zwei freie Parameter, die wir im folgenden anpassen wollen:
Mit Hilfe der Shooting Methode kann man Randwertprobleme numerisch auf Anfangswertprobleme zurückführen. Es wird zu Anfang ein Intervall geraten in dem die mögliche Lösung liegt. Für die Mitte des Intervalls wird das Problem numerisch gelöst und mit dem Randwert verglichen. Es wird nun mit Hilfe des Bisektionverfahrens das Intervall solange halbiert bis die Lösung ausreichend genau gefunden wurde.
Wir modifizieren die Massen der Flutberge so, dass sie für die aktuellen Zunahme der Tageslänge pro 100 Jahre mit dem oben angegebenen Literaturwert $\tau$ übereinstimmt (hier bleibt $k$ unverändert).
In der Differentialgleichung für die Erdrotationsbeschleunigung sehen wir, dass diese linear mit der Masse der Flutberge zusammenhängt. Die berechnete Erdrotationsbeschleunigung ist zu schnell, deshalb müssen wir die Masse der Flutberge verkleinern. Wir wählen deshalb für die Shooting Methode als Suchintervall Null bis zur aktuellen Ozeanmasse. Es wäre möglich die Massen der Flutberge einzeln zu variieren, wir haben uns aber entschieden sie für die folgenden Berechnungen gleich zu setzen.
Die verwendete Ozeanmasse entspricht bereits den reelen Begebenheiten, deshalb ist es nicht sinnvoll diese zu fitten. Alternativ können wir die Reibungskonstante verändern. Auch sie ist aus den gleichen Gründen zu groß und wir wählen als Suchintervall Null bis zum aktuellen Wert. Die Masse der Flutberge lassen wir hier unverändert.
Zum Ausführen der Shooting-Methode müssen wir noch eine Toleranz angeben. Der Literaturwert $\tau = 0.0021$ ist mit einer Genauigkeit von $10^{-5}$ angegeben, deshalb wählen wir diese als Toleranz.
def rotation_perdiod_change_mOcean(mOcean: float):
'''Calculates the change in rotation_period after 100 years, given a value for mOcean'''
initial_values = iv_stable_orbit_4body_friction()
# Change the mass of the two tides to each be half the value of mOcean for the simulation
initial_values[-1][-1], initial_values[-1][-2] = mOcean/2, mOcean/2
_, _, _, _, _, _, _, _, _, _, omega_E = four_body_problem_friction(*initial_values, k=2*10**(-12), t_max=100*365*24*3600) # Simulate 100 years
return (2*pi/omega_E[-1]) - (2*pi/omega_E[0]) # Calculate the change in rotation period as described above
def rotation_period_change_k(k: float):
'''Calculates the change in rotation_period after 100 years, given a value for k'''
_, _, _, _, _, _, _, _, _, _, omega_E = four_body_problem_friction(*iv_stable_orbit_4body_friction(), k=k, t_max=100*365*24*3600) # Simulate 100 years
return (2*pi/omega_E[-1]) - (2*pi/omega_E[0]) # Calculate the change in rotation period as described above
def shooting_method(func, target: float, val_min: float, val_max: float, tol: float):
'''Find the value of x for which func(x) is in the interval [target-tol, target+tol]
Search in the interval [val_min, val_max]'''
# Check if the solution is in the passed interval
if func(val_min) > target or func(val_max) < target:
raise ValueError('The target is not in the passed interval')
print("Target in interval")
# Iterate until the solution is within the tolerance
while True:
x = (val_min + val_max)/2 # Calculate the midpoint of the interval
func_val = func(x) # Calculate the value of the function at the midpoint
if np.abs(func_val- target) < tol: # Solution found if the absolute distance to the solution is smaller than the tolerance
return x
elif func_val < target: # The lower bound can be raised
val_min = x
else: # The upper bound can be lowered
val_max = x
# Fit the value of mOceans so that the change in rotation period is tau
if not fastExecution:
mOzean_fit = shooting_method(func=rotation_perdiod_change_mOcean, target=tau, val_min=0, val_max=mOzean, tol=0.00001)
else:
mOzean_fit = 2.13623046875e+16
# Fit the value of k so that the change in rotation period is tau
if not fastExecution:
k_fit = shooting_method(func=rotation_period_change_k, target=tau, val_min=0, val_max=2*10**(-12), tol=0.00001)
else:
k_fit = 3.06640625e-17
print(f'mOzean_fit: {mOzean_fit:.3g} kg')
print(f'k_fit: {k_fit:.3g} 1/m')
mOzean_fit: 2.14e+16 kg k_fit: 3.07e-17 1/m
Wir erhalten wie erwartet durch den Fit für Ozeanmasse und Reibungskonstante kleinere Werte als unsere Startwerte.
Der gefittete Wert für die Ozeanmasse liegt aber mit $2.14 \cdot 10^{16} \mathrm{kg}$ deutlich unter dem vorher verwendeteten Literaturwert von $ 1.4 \cdot 10^{21} \mathrm{kg} $ Quelle: Ocean#Weigth.
Die Abweichungen lassen sich durch mehrere Näherungen in der Simulation erklären:
Nach den obigen Betrachtungen scheint es also sinnvoller statt der Ozeanmasse die Reibungskonstante zu fitten. Diese verhält sich auch viel mehr wie ein freier Parameter, da sie viele in der Realität auftretende Effekte zusammenfasst. (z.B. Landmassen, verschiedene Meerestiefen, Strömungen)
if not fastExecution:
# 4 body with friction for 100 years: 11min
t_kfit, x_E_kfit, y_E_kfit, x_M_kfit, y_M_kfit, phi_1_kfit, phi_2_kfit, phi_E_kfit, omega_1_kfit, omega_2_kfit, omega_E_kfit = four_body_problem_friction(*iv_stable_orbit_4body_friction(), k=k_fit, t_max=100*365*24*3600, rtol=1e-12, atol=1e-12)
x_1_kfit, y_1_kfit, x_2_kfit, y_2_kfit = x_E_kfit + RErde*np.cos(phi_1_kfit), y_E_kfit + RErde*np.sin(phi_1_kfit), x_E_kfit + RErde*np.cos(phi_2_kfit), y_E_kfit + RErde*np.sin(phi_2_kfit)
if not fastExecution:
fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(15, 10))
plt.rc('axes', axisbelow=True)
plt.suptitle("Verlgeich 4 Körper Problem $k=2*10^{-12}$ / $k=3*10^{-17}$")
def plot_subplot(ax, lines, title, xlabel, ylabel):
for x, y, label, color in lines:
ax.plot(x, y, color, linewidth=1, label=label)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_title(title)
ax.grid('minor', 'minor', linestyle='-', linewidth=0.2)
ax.grid('major', 'major', linestyle='-', linewidth=0.8)
ax.minorticks_on()
ax.legend()
# Distance from the origin for earth and moon
r_Earth_fl = np.sqrt(x_E_fl**2 + y_E_fl**2)
r_Moon_fl = np.sqrt(x_M_fl**2 + y_M_fl**2)
r_Earth_kfit = np.sqrt(x_E_kfit**2 + y_E_kfit**2)
r_Moon_kfit = np.sqrt(x_M_kfit**2 + y_M_kfit**2)
plot_subplot(axes[0][0], [(x_E_fl*10, y_E_fl*10, 'Erde (*10)', 'k'), (x_M_fl, y_M_fl, 'Mond', 'r'), (x_1_fl*10, y_1_fl*10, 'Flutberg 1 (*10)', 'b'), (x_2_fl*10, y_2_fl*10, 'Flutberg 2 (*10)', 'b')], "Bahnen des 4-Körper Problems $k=2*10^{-12}$", "x / m", "y / m")
plot_subplot(axes[0][1], [(t_fl, r_Moon_fl, 'Mond', 'r')], "Abstand des Mondes zum Ursprung $k=2*10^{-12}$", "t / s", "r / m")
plot_subplot(axes[0][2], [(t_fl, r_Earth_fl, 'Erde', 'k')], "Abstand der Erde zum Ursprung $k=2*10^{-12}$", "t / s", "r / m")
plot_subplot(axes[1][0], [(x_E_kfit*10, y_E_kfit*10, 'Erde (*10)', 'k'), (x_M_kfit, y_M_kfit, 'Mond', 'r'), (x_1_kfit*10, y_1_kfit*10, 'Flutberg 1 (*10)', 'b'), (x_2_kfit*10, y_2_kfit*10, 'Flutberg 2 (*10)', 'b')], "Bahnen des 4-Körper Problems $k=3*10^{-17}$", "x / m", "y / m")
plot_subplot(axes[1][1], [(t_kfit, r_Moon_kfit, 'Mond', 'r')], "Abstand des Mondes zum Ursprung $k=3*10^{-17}$", "t / s", "r / m")
plot_subplot(axes[1][2], [(t_kfit, r_Earth_kfit, 'Erde', 'k')], "Abstand der Erde zum Ursprung $k=3*10^{-17}$", "t / s", "r / m")
plt.tight_layout()
plt.savefig("4body_friction_comp1_k_fit.png", dpi=300)
plt.close()
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 5))
plt.rc('axes', axisbelow=True)
plt.suptitle("Verlgeich 4 Körper Problem $k=2*10^{-12}$ / $k=3*10^{-17}$")
plot_subplot(axes[0][0], [(t_fl, omega_1_fl, '$\omega_1$ Flutberg 1', 'k'), (t_fl, omega_2_fl, '$\omega_2$ Flutberg 2', 'r'), (t_fl, [2*pi/TMondBahn]*np.shape(t_fl)[0], 'Startwinkelgeschwindigkeit Flutberg', 'g')], "Winkelgeschwindigkeiten der Flutberge $k=2*10^{-12}$", "t / s", "$\omega$ / rad/s")
plot_subplot(axes[0][1], [(t_fl, omega_E_fl, '$\omega_E$ Winkelgeschwindigkeit Erde', 'k'), (t_fl, [2*pi/TErdRotation]*np.shape(t_fl)[0], 'Startwinkelgeschwindigkeit Erde', 'g')], "Winkelgeschwindigkeit der Erde $k=2*10^{-12}$", "t / s", "$\omega$ / rad/s")
plot_subplot(axes[1][0], [(t_kfit, omega_1_kfit, '$\omega_1$ Flutberg 1', 'k'), (t_kfit, omega_2_kfit, '$\omega_2$ Flutberg 2', 'r'), (t_kfit, [2*pi/TMondBahn]*np.shape(t_kfit)[0], 'Startwinkelgeschwindigkeit Flutberg', 'g')], "Winkelgeschwindigkeiten der Flutberge $k=3*10^{-17}$", "t / s", "$\omega$ / rad/s")
plot_subplot(axes[1][1], [(t_kfit, omega_E_kfit, '$\omega_E$ Winkelgeschwindigkeit Erde', 'k'), (t_kfit, [2*pi/TErdRotation]*np.shape(t_kfit)[0], 'Startwinkelgeschwindigkeit Erde', 'g')], "Winkelgeschwindigkeit der Erde $k=3*10^{-17}$", "t / s", "$\omega$ / rad/s")
plt.tight_layout()
plt.savefig("4body_friction_comp2_k_fit.png", dpi=300)
plt.close()
# Display image set the size to be 500x500 pixels
display(Image(filename='4body_friction_comp1_k_fit.png', width=1000, height=500))
display(Image(filename='4body_friction_comp2_k_fit.png', width=1000, height=500))
Da die gefittete Reibungskonstante um fünf Größenordnungen kleiner ist, sind die zuvor diskutierten Reibungseffekte sehr viel kleiner. Die Tageslängenänderung lässt sich jedoch trotzdem erkennen. Die anderen Effekte (Entfernung des Mondes von der Erde, Einschwingen), sind jedoch so klein, dass sie in den Plots nicht zu sehen sind.
Nun betrachten wir das Erde-Mond System mit $N$ Massepunkten auf der Erdoberfläche, welche die Wassermassen darstellen. Ziel ist es das Ausbilden der Flutberge zu simulieren. Für die Mond- und Erdkoordinaten verwenden wir die Differentialgleichung des Zweikörperproblems, als Näherung um den Rechenaufwand zu reduzieren. Die Näherung wird später genauer diskutiert. Die Kraft auf die Masspunkte ist gleich zur Kraft auf die zwei Flutberge mit Reibung, welche oben erläutert wurde. Die Differentialgleichung für $\omega_{\mathrm E}$ wird mit $N$ Massepunkten, analog zum Fall mit zwei Flutbergen, zu: $\newline$ \begin{equation*} \dot{\omega}_{\mathrm E} = \frac{5kR_{\mathrm E}}{2m_{\mathrm E}} \sum_{i = 1}^{N} m_i\lvert \omega_i - \omega_{\mathrm E} \rvert(\omega_i - \omega_{\mathrm E}) \end{equation*}
N = 25 # Amount of tides
Wir verwenden wenn möglich die bereits implementierten Gleichungen. Die Winkel und Winkelbeschleunigungen der Flutberge werden als Liste übergeben, um $N$ leicht änden zu können.
def eq_motion_Nbody(t, state, mass, k):
'''This function calculates the derivatives of the state vector the 2 body problem with N Tides to be passed to solve_ivp. \\
state: state vector has the following structure: [x_E, y_E, x_M, y_M, phi_E, vx_E, vy_E, vx_M, vy_M, omega_E, phi_1, phi_2, ..., phi_N, omega_1, omega_2, ..., omega_N]
mass: list of masses e.g. [m_E, m_M, m1, m2, ..., mN]'''
# Unpack the state vector and masses:
x_E, y_E, x_M, y_M, phi_E, vx_E, vy_E, vx_M, vy_M, omega_E, phi_tide, omega_tide = *state[0:10], [state[10+i] for i in range(N)], [state[10+N+i] for i in range(N)]
m_E, m_M, m_tide = *mass[0:2], [mass[2+i] for i in range(N)]
# Earth and moon acceleration from the 2 body problem:
_, _, _, _, ax_E, ay_E, ax_M, ay_M = eq_motion_2body(t, [x_E, y_E, x_M, y_M, vx_E, vy_E, vx_M, vy_M], [m_E, m_M])
# Angular momentum of the earth:
alpha_E = (5*k*RErde)/(2*m_E) * np.sum([m_tide[i]*abs(omega_tide[i] - omega_E)*(omega_tide[i] - omega_E) for i in range(N)])
# The tides' cartesian coordinates:
x_tide, y_tide = [x_E + RErde*np.cos(phi_tide[i]) for i in range(N)], [y_E + RErde*np.sin(phi_tide[i]) for i in range(N)]
# The earth's angular acceleration:
alpha_tide = [tides_acceleration_friction(x_M, y_M, x_tide[i], y_tide[i], phi_tide[i], ax_E, ay_E, omega_tide[i], omega_E, m_M, k) for i in range(N)]
return [vx_E, vy_E, vx_M, vy_M, omega_E, ax_E, ay_E, ax_M, ay_M, alpha_E] + omega_tide + alpha_tide
Da für die Erd- und Mondkoordinaten die Zweikörper-Differentialgleichug verwendet wird wählen wir für diese dieselben Anfangsbedingungen. Die Massepunkte starten gleichmäßig verteilt auf der Erdoberfläche mit der mittleren Winkelgeschwindigkeit der Flutberge mit Reibung.
if not fastExecution:
_, _, _, _, _, _, _, _, omega_Niv, _, _ = four_body_problem_friction(*iv_stable_orbit_4body_friction(), k=k_fit, t_max=100*365*24*3600, rtol=1e-12, atol=1e-12)
omega_i0 = np.mean(omega_Niv)
else:
omega_i0 = 2.7107783823217206e-06
def iv_stable_orbit_Nbody():
'''Initial conditions for the Earth-Moon-NTides system in a stable orbit'''
# Each of the tides gets an nth of the oceans mass
mass = [mErde, mMond] + [mOzean/N for i in range(N)]
# Because the DEQ for the 2 body problem is used for the earth and moon movement, we can reuse the initial conditions (Tide mass evenly distributed around the earth)
x0_Erde, x0_Mond, v0_Erde, v0_Mond, _ = iv_stable_orbit_2body()
_, _, _, phi0_E, _, _, _, omega0_E, _ = iv_stable_orbit_4body_friction()
phi0 = [(2*pi*i)/N for i in range(N)] # The tides all start spaced evenly around the earth
#omega0 = [2*pi/TMondBahn]*N # The tides all start with the same angular velocity as the moon
omega0 = [omega_i0]*N
return [x0_Erde, x0_Mond, phi0_E, v0_Erde, v0_Mond, omega0_E, phi0, omega0, mass]
def N_body_problem(pos_body_1: list, pos_body_2: list, phi_E: float, vel_body_1: list, vel_body_2: list, omega_E: float, phi: list, omega: list, mass: list, k: float, t_max: float, t_start=0):
solution = solve_ivp(fun=eq_motion_Nbody, t_span=[t_start, t_max], y0=[*pos_body_1, *pos_body_2, phi_E, *vel_body_1, *vel_body_2, omega_E, *phi, *omega], args=(mass, k), rtol=1e-12, atol=1e-12)
x_E, y_E, x_M, y_M, phi_E, vx_E, vy_E, vx_M, vy_M, omega_E = solution.y[0:10]
phi_tide = [solution.y[10+i] for i in range(N)]
omega_tide = [solution.y[10+N+i] for i in range(N)]
return [solution.t, x_E, y_E, x_M, y_M, phi_tide, omega_E, omega_tide]
t, x_E, y_E, x_M, y_M, phi_tide, omega_E, omega_tide = N_body_problem(*iv_stable_orbit_Nbody(), k=k_fit, t_max=20*TMondBahn)
# Cartesian coordinates of the tides:
x_tide, y_tide = [x_E + RErde*np.cos(phi_tide[i]) for i in range(N)], [y_E + RErde*np.sin(phi_tide[i]) for i in range(N)]
# Plot the movement of the earth, moon and tides
plt.rc('axes', axisbelow=True)
plt.figure(figsize=(5, 5))
plt.title("Bahnen des N-Flutberg-System")
plt.grid('minor', 'minor', linestyle='-', linewidth=0.2)
plt.grid('major', 'major', linestyle='-', linewidth=0.8)
plt.minorticks_on()
plt.xlabel("x")
plt.ylabel("y")
plt.plot(x_M/10, y_M/10, 'r', linewidth=1, label='Mond (/10)')
for i in range(N):
plt.plot(x_tide[i], y_tide[i], 'b', linewidth=0.5, alpha=0.2)
plt.plot(x_E, y_E, 'k', linewidth=1, label='Erde')
plt.legend()
plt.show()
plt.close()
# Plot the angles of the tides over time, subtracting the angle of the moon
plt.rc('axes', axisbelow=True)
plt.figure(figsize=(5, 5))
plt.title("Winkel der Flutberge über Zeit, abzüglich $t\cdot \omega_0$")
plt.grid('minor', 'minor', linestyle='-', linewidth=0.2)
plt.grid('major', 'major', linestyle='-', linewidth=0.8)
plt.minorticks_on()
plt.xlabel("$t$ / s")
plt.ylabel("$\phi - t\cdot \omega_0$ / rad")
for i in range(N):
plt.plot(t, phi_tide[i] - t*2*pi/TMondBahn, 'k', linewidth=0.7, alpha=0.7)
plt.show()
plt.close()
# Plot the angular velocity of the earth over time
plt.rc('axes', axisbelow=True)
plt.figure(figsize=(5, 5))
plt.title("Winkelgeschwindigkeit der Erde über Zeit")
plt.grid('minor', 'minor', linestyle='-', linewidth=0.2)
plt.grid('major', 'major', linestyle='-', linewidth=0.8)
plt.minorticks_on()
plt.xlabel("$t$ / s")
plt.ylabel("$\omega_E$ / rad")
plt.plot(t, omega_E, 'k', linewidth=1, label='$\omega_E$')
plt.show()
plt.close()
# Animation of the N body problem
if not fastExecution:
fig = plt.figure(figsize=(5, 5))
ax = plt.axes(xlim=(-4e8, 4e8), ylim=(-4e8, 4e8))
ax.set_title("Animation des N-Flutberg-Systems")
earth_line, = ax.plot([], [], marker='o', lw=0.1, color='blue', label='Erde (*10)')
moon_line, = ax.plot([], [], marker='o', lw=0.1, color='red', label='Mond')
tide_lines = [ax.plot([], [], marker='o', lw=0.1, color='black', label='Flutberg {i} (*10)', alpha=0.15)[0] for i in range(N)]
def animate(i):
earth_line.set_data([x_E[i]*10], [y_E[i]*10])
moon_line.set_data([x_M[i]], [y_M[i]])
for j in range(N):
tide_lines[j].set_data([x_tide[j][i]*10], [y_tide[j][i]*10])
return earth_line, moon_line, *tide_lines
anim = FuncAnimation(fig, animate, init_func=None, frames=t.shape[0], interval=30, blit=True)
anim.save('Erde_Mond_NFlutberge.gif', writer='pillow')
plt.close()
display(Image(data=open('Erde_Mond_NFlutberge.gif','rb').read(), format='png'))